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An iterative algorithm for the reconstruction of an 
unknown quantum state from the results of incompati- 
ble measurements is proposed. It consists of Expectation- 
Maximization step followed by a unitary transformation of 
the eigenbasis of the density matrix. The procedure has been 
applied to the reconstruction of the entangled pair of photons. 



I. INTRODUCTION 

Predictions of advanced theories are more and more 
complex and more and more accurate. This may be rec- 
ognized in recent progress of quantum theory. In many 
apphcations there is a need to determine the quantum 
state of the system. For this purpose the quantum to- 
mography has been developed. There is an extended bib- 
hography concerning this topic covering various fields of 
possible applications y. However, all the varied and pre- 
cise experiments which have been carried out over many 
years and which rely on quantum physics have not need 
the total amount of information coded in quantum state. 
This fact is reflected in the standard treatment of quan- 
tum tomography. When the standard quantum tomog- 
raphy is adopted for reconstruction of the state from re- 
alistic noisy data, one often runs into unphysical results. 
Standard methods are also known to be prone to cre- 
ation of various artifacts in the reconstructed state. All 
these flaws are usually paid little or no attention in the 
scientific literature. They are simply being regarded as 
unavoidable errors of reconstructions, which fall within 
the corresponding "error bars" . Here we would like to 
stress that the mentioned drawbacks of the standard to- 
mographic techniques are actually much more serious, es- 
pecially if the reconstructed state is to be of further use. 
There might be no clue as to which modifications to the 
reconstructed "state" should be done in order to make 
its density operator semi-positive definite and retain the 
efficiency of the reconstruction as high as possible at the 
same time. This seems to be crucial for the potential 
application in quantum information. To quantify a frag- 
ile effect of entanglement, various entropic principles are 
used 1^. This feature is sensitive to semi-positive defi- 
nitcness of the reconstructed state - a necessary condition 
of successful reconstruction. 

The purpose of this Rapid Communication is twofold. 
At first a simple iterative algorithm for maximum like- 
lihood (ML) estimation of quantum state resembling 
"climbing the hill of the likelihood" will be derived. This 
result may be easily implemented numerically, and inter- 
preted in quantum theory as generalized measurement. 
The algorithm will be illustrated on the example of re- 



construction of entangled state, representing an impor- 
tant example in quantum information processing. 

Let us illustrate the motivation of quantum tomog- 
raphy considering the repeated measurement. Assume 
that we are given a finite number N of identical samples 
of the system, each in the same but unknown quantum 
state described by the density operator p. Given those 
systems our task is to identify the unknown true state 
p from the results of measurements performed on them 
as accurate as possible. For simplicity we will assume 
sharp measurements in the sense of von Neumann. As a 
result of each measurement the state of the input system 
is projected into a pure state, which is the reading of the 
measuring apparatus. Let us assume, for concreteness, 
that M different outcomes of measurements have been 
observed. The relative frequencies fj of occurrences of 
the observed results 
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then comprise the data which the true state p should be 
inferred from. For the sake of simplicity, the measure- 
ment performed will be assumed as complete, i.e. 



H^J2i\y,){y,\}^l. 



This condition will be released later to the case of incom- 
plete measurements. 

The probabilities of occurrences of various outcomes 
are generated by the true quantum state p according to 
the well-known handy quantum rule 



Pj = {yMvj)- 



(2) 



If the probabilities pj of getting a sufficient number of 
different outcomes \yj) were known, it would be possible 
to determine the true state p directly by inverting the 
linear relation (ph. This the philosophy behind the "stan- 
dard" quantum tomographic techniques. For example, in 
the rather trivial case of a spin half particle, the prob- 
abilities of getting three linearly independent projectors 
determine the unknown state uniquely. Here, however, a 
serious problem arises. Since only a finite number of sys- 
tems can be investigated, there is no way how to find out 
these probabilities. The only data one has at disposal are 
the relative frequencies fj, which sample the principally 
unknowable probabilities pj . It is obvious that for small 
number of runs the true probabilities pj and the corre- 
sponding detected frequencies fj may differ substantially. 
As a result of this, the modified realistic problem 



fj = ivMyj) 



(3) 



has generally no solution on the space of semi-positive 
definite hermitian operators describing physical states. 

Probabilistic interpretation of quantum theory sug- 
gests that a sort of statistical treatment of the observed 
data might be more natural and appropriate than the de- 
terministic treatment described above. The philosophy 
of the reconstruction method differs from the philosophy 
of standard methods. The basic question of the stan- 
dard methods: "What quantum state is determined by 
the measured data?" is replaced with a more modest 
one: "What quantum state is most likely in view of the 
measured data?"; this seems to be in accordance with 
the probabilistic interpretation of the quantum theory 
Q . More specifically, instead of trying to invert the lin- 
ear relation (j^, we look for a density operator pe, which 
generates through Eq. (0) probabilities pj that are as 
"close" to the observed frequencies fj as possible, i.e. 



Pe = arg{mind[f,p(p)]}, 
p 



(4) 



where d[., .] stands for some measure of distance between 
the two probability distributions p and data f . At first 
sight it might seem that there is no reason to prefer one 
particular metric to another one - different metrics lead- 
ing to different results through Eq. (0). This ambiguity 
can be resolved by considering the formal description of 
the reconstruction process M . If the whole measurement 
and subsequent reconstruction is looked at as a single 
generalized measurement, then the relation between the 
actually performed measurement and resulting probabil- 
ity operator measure becomes particularly simple and 
easy to interpret for the metric known as the relative 
entropy or KuUback-Leibler divergence ||] : 



4^: P]^-Y1 f'J l^Pj 



(5) 



Solving the problem (0) with the metric (ra) is equivalent 
to finding the maxima of the likelihood functional 



t^ip) = \{{yMvj) 



f: 



(6) 



Thus we are led to the Maximum Likelihood principle as 
the preferred way of doing the quantum state reconstruc- 
tion. 

ML methods are well-known in the field of inverse 
problems and they have found many applications in re- 
constructions and estimations so far [|6|. Unfortunately, 
except in most simple cases, the maximization of the like- 
lihood functional is a challenging problem on its own. 
A necessary condition for an extreme of the likelihood 
functional (Jg) can be derived in the form of the nonlin- 
ear operator equation for the density matrix p |^,0 and 
this equation may be interpreted as closure relation for a 
quantum measurement. In the classical signal processing 
an important role is played by the so called linear and 
positive (LinPos) problems ||] . Since these are closely re- 
lated to the problem of quantum state reconstruction it is 



worthwhile to recall how the positive and linear problems 
can be dealt with using the ML approach. 

Let us consider that the probabilities pj of getting out- 
comes yj are given by the following linear and positive 
relation 



Pj 



/ , '"'i'T-ij^ 



p,r, h > 0. 



(7) 



Here r is the vector describing the "state" of the system. 
For example, the reconstruction of a one-dimensional 
object from the noiseless detection of its blurred im- 
age could be accomplished by inverting the relation (|^), 
where r and p would be the normalized intensities of 
the object and image, and h would describe the blurring 
mechanism. Again here the presence of noise {fj ^ pj) 
tends to spoil the positivity of the reconstructed intensity 
r. 

The solution to LinPos problems in the sense of (H) 
can be found using the Expectation-Maximization (EM) 
algorithm Jpj. In the case of the discrete one-dimensional 
problem (^ the unknown object r is reconstructed by 
means of the following iterative algorithm S : 
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^p^.(l.(«-l))^ 
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which is initialized with a positive vector r (r^ > OVz). 

The iterative algorithm (pi) for solving the LinPos prob- 
lems is convenient from the point of view of the numerical 
analysis. It is certainly much more convenient than the 
direct multidimensional maximization of the correspond- 
ing ML functional InC — J2j fj ^"^Pj @- This brings us 
back to the problem of quantum state reconstruction. It 
would be nice to have a similar iterative algorithm for 



dealing with the problem 
imizing the ML functional 



[or equivalently for max- 
t)]. On the one hand it is 



clear that the problem of quantum state reconstruction 
is not a linear and positive problem, since the quantum 
rule (H) cannot be rewritten to the form of Eq. (M) with a 
known positive kernel h. As a consequence of this the EM 
algorithm cannot be straightforwardly applied here. On 
the other hand the reconstruction of the elements of the 
density matrix becomes a LinPos problem if the eigen- 
basis diagonalizing the density matrix is known. In this 
case the unknown density matrix can be parametrized as 
follows 



P = ^rk\(l)k){<l>k 



p\4>k) = rk\(t>k) 



(9) 



where r^ are eigenvalues of p, the only parameters which 
remain to be determined from the performed measure- 
ment. Using the parameterization (M) the quantum rule 
(g) may be easily rewritten to the form of LinPos problem 
Eq. (0). ^ 

This hints on splitting the quantum state reconstruc- 
tion into two subsequent steps: the reconstruction of the 
eigenvectors of p in a fixed basis, which represents the 



classical part of the problem, followed by the "rotation" 
of the basis {|(?!>z)} in the "right" direction using the uni- 
tary transformation 



W,)m = u\<Pk){Mu^. 



Its infinitesimal form reads 



U = e''^ 



1 



eG. 



(10) 



(11) 



Here G is a Hermitian generator of the unitary trans- 
formation and e is a positive real number which is small 
enough in order to make the second equality in (O) ap- 
proximately satisfied. 

Consider now the total change of log- likelihood caused 
by the change of diagonal elements of density matrix and 
rotation of basis. Keeping the normalization condition 
Trp — 1, the first order contribution to the variation 
reads 



k 

+zeTr{G[p,i?]}. 



(12) 



The operator R appearing here plays an important role 
in this treatment. It is semi-positively definite Hermitian 
operator comprising results of the measurement 



R 



3 ■' 



(13) 



Notice this operator depends on the old density matrix 
p through Eq. (|). 

Inspection of Eq. (|lj) reveals a simple strategy how to 
make the likelihood of the new state p' as high as possible 
[within limits of the validity of the linearization (|ll|), of 
course]. In the first step the first term on the right hand 
side of Eq. (jlj) is maximized by estimating the eigen- 
values of the density matrix keeping its eigenvectors \(j)k) 
constant. The iterative algorithm (g) described above 
can straightforwardly be applied to this LinPos prob- 
lem. In the second step the likelihood can further be 
increased by making the second term on the right hand 
side of Eq. ( p^ positive. This is accomplished by a suit- 
able choice of the generator of the unitary transformation 
( TWj . Reminding the norm induced by the scalar product 
defined on the space of operators , {A,B) = TrjA^i?}, 
the generator G may be chosen as 



G^i[p,R]. 



(14) 



Its form guarantees the non-negativity of the contribu- 
tion to the likelihood. This choice is also optimal in the 
sense of the above introduced scalar product. 

Now we have at our disposal all ingredients comprising 
the EMU quantum state reconstruction algorithm which 
represents the main result of the present article. Start- 
ing from some positive initial density matrix p this esti- 
mate is improved, first by finding new eigenvalues using 



the EM iterative algorithm (g), and then again by find- 
ing new eigenvectors by unitary (U) transforming the old 
ones according to Eqs. (p^pl]) and (|l^. These two steps 
are repeated. Continued repetition of the two steps, each 
monotonically increasing the likelihood of the current es- 
timate, resembles of climbing a hill. 

The proposed EMU algorithm naturally leads to the 
previously introduced extremal equation for the density 
matrix am . The stationary point of EMU algorithm is 
characterized by the vanishing variance of the log likeli- 
hood (O) . Since the variations Srk , e are arbitrary pa- 
rameters, this is equivalent to the Lagrange-Euler equa- 
tion for density matrix 



Rpe 



(15) 



This nonlinear operator equation has recently been de- 
rived using the variational principle in Ref. 0| and using 
the inequalities Q. The EMU algorithm presented here 
thus provides us with a different route to its derivation, 
which is perhaps more appealing from the physical point 
of view, and suitable for implementation of numerical al- 
gorithm. Notice, however, that this "parameter estima- 
tion" may be interpreted as a generalized measurement 
since i? := 1 on the space where the reconstruction has 
been done. 

These results should be modified in the case of incom- 
plete detection. Provided that H ^ 1, the closure rela- 
tion may be always recovered in the form 



J2H''/my,}{y,\}H-'/'^l. 



(16) 



This corresponds to the renormalization of the probabili- 
ties pj = {yj\p\yj) in the likelihood (o) to the normalized 
probabilities pj -^ Pjl^iPi- This formulation incor- 
porates the case of incomplete detection. Notice, that 
extremal equation possesses again the form of equation 
( |l5| ) for the renormalized quantities 



R^R' = {H')-^I^R{H')-^'^, 
p,-.p',^{H'f''p,{H'f'\ 



E,Pj 



(17) 



All the conclusions derived for complete measurements 
may be extended to this case of incomplete measurement 
as well. This formulation coincides with the estimation 
provided that assumption of Poissonian statistics is used. 
Assume that rii samples the mean number of particles 
npi, where pi is as before the prediction of quantum the- 
ory for detection the i-th channel and n is unknown mean 
number of particles. The relevant part of log-likelihood 
corresponding to the Poissonian statistics reads 



In £ oc 2_, "-i lii("-Pi) — n /^ Pi 



(18) 



The extremal equation for n may be easily formulated as 
the condition n = J2i''^i/J2iPi- Inserting this estimate 
of unknown mean number of Poissonian particles into 
the log-likelihood reproduces the renormalized likelihood 
function. 

The proposed EMU algorithm has been applied to 
the reconstruction of two-photon entangled state gener- 
ated by the spontaneous downconversion source of White 
et al. flQ]. White et al. measured the nominal state 
\HH) + \VV) along sixteen distinct directions: {\yj)} = 
{\HH), \HV), \VH), \VV), \HD), \HL), \DH), \RH), 
\DD), \RD), \RL), \DR), \DV), \RV), \VD), \VL)}. 
H, V, D, R, and L being horizontal, vertical, diagonal, 
right circular, and left circular polarization, respectively. 
Counted numbers of coincidences along these directions 
can be found in pvj \. 

We have used the experimental data together with the 
proposed algorithm to reconstruct the true state gener- 
ated by the source of entangled photon pairs. Due to 
various sources of errors the true state is expected to dif- 
fer from the nominal state. Notice that the chosen mea- 
surements are not complete, that is ^ \yj){yj\ does not 
represents the resolution of unity. This has been taken 
into account, see Eq. (p7|). 

Starting from maximally mixed state {\HH){HH\ + 
\VV){VV\ + \HV){HV\ + \VH){VH\)/'i,, new eigenval- 
ues and eigenvectors of density matrix are found using 
Eqs. (H) and ([lO|). This has been repeated until a sta- 
tionary point of the iteration process has been attained. 
The diagonal representation of the reconstructed density 
matrix reads 



pf ^ = 0.962 |0i}((/)i| + 0.038 \c^2){h 



(19) 



The other two eigenvalues are zero. The eigenvectors |0i) 
and 102) are given in Tab. (|). 

The reconstructed density matrix (|l9|) agrees well with 
the qualitative reasoning given in [|10[. Namely, the re- 
constructed state is almost pure state - slightly rotated 
nominal state \HH) + \VV). The apparent incompati- 
bility of the nominal state with the registered data was 
interpreted in [nO| as the result of possible slight misalign- 
ments of the axes of analysis systems with respect to the 
axes of the downconversion source. This is, of course, re- 
flected in the reconstructed state ([19|), which quantifies 
such misalignments and might serve for hunting down the 
errors, and calibrating the experimental setup. 

Notice also that the reconstructed density matrix (|l9| ) 
is semi-positive definite. This should be contrasted with 
the result of standard reconstruction. Direct inversion 
of Eq. (0) yields the density matrix having the following 
diagonal representation [|ll| 

ri = 1.022, 7-2 ^ 0.068, rg = -0.065, r^ = -0.024 (20) 

The corresponding eigenvectors need not be specified 
here. Apparently, direct inversion (standard tomogra- 
phy) leads to unphysical non-positive definite result. It 



is worth to notice that the negative eigenvalues are com- 
parable in magnitude with non-diagonal elements of p^^ 
in H-V basis, see Tab. fl This is a nice example of sit- 
uation when standard methods fail even though rather 
high number of particles (tens of thousands) has been 
registered. ML reconstruction provides always physically 
sound results. Moreover, it represents genuine quantum 
measurement of entangled state. 
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\VV) 
\VH) 
\HV) 
\HH) 



0.696 - 0.027i 
-0.050 - 0.020i 
-0.040 + 0.015i 

0.712 - 0.062J 



0.630 + 0.071J 
-0.284 -F0.174J 
-0.150 -0.247J 
-0.634 - 0.035J 



TABLE L Eigenvectors of the reconstructed density matrix. 



